In many areas of digital signal processing, there is a need to generate a sequence of sampled sine and cosine functions. Examples that come immediately to mind are Fourier transform evaluations, digital quadrature demodulators, and generation of test tones.
A direct method is to use a stored table look-up procedure wherein a table of function values corresponding to argument values is provided. This method offers high speed, and as much accuracy as desired. Its main problem is excessive memory storage. If a 1024 point sequence is required for example, with 16 bit arithmetic, a 32 Kbyte ROM is needed just to store the sine and cosine waveforms (assuming no repeats of the data).
In a second method, if we are willing to invest some computation time per cycle, this memory requirement can be reduced dramatically by providing an arithmetic unit for calculating function values. For example, the simple recursion EQU sin [(n+1)x]=cos (x) sin (nx)+sin (x) cos (nx) EQU cos [(n+1)x]=cos (x) cos (nx)-sin (x) sin (nx)
requires four multiplications per cycle and needs only a few memory locations.
One has an uneasy feeling about using these formulas however, as an extrapolation is involved. Round-off errors will occur, and grow as the extrapolation proceeds. The related phasor tends to spiral inwards or outwards, rather than retain constant amplitude.
One known way of overcoming this extraplation problem is to slightly lengthen the phasor increment. The running phasor is thus given a small outward movement on each update. Then, assuming the arithmetic used in the dsp calculating device saturates for sample values greater than 1, the output will be (slightly) clipped, but stable. The problem with this method is that it can build up a considerable error before clipping occurs, so the resulting error is large. An alternative, third method is based on interpolation (1), (2). The sine and cosine functions, initially accurately calculated over a small set of sample points, which span the range of interest, are stored. These values are used as inputs to an algorithm performed by an arithmetic unit, which uses interpolation calculations to find all remaining values. Small errors can and will occur, but they are significantly smaller than errors from an extrapolation algorithm.
The price paid for this increased accuracy is a small increase in computation time over the first method, and a modest increase in memory size over the second.
In one simple method of interpolation, described in (3), a short table of length (log.sub.2 N)+2 where N is the number of derived angular values (spaced at 2.pi./N) including sec values of (.pi./4), (.pi./8), (.pi./16) . . . (.pi./N), is used to calculate successive levels of a table of interpolated function values, starting from values at .pi./4, 0 and interpolating to find values at .pi./8 and 3.pi./8, then interpolating between these to find those at .pi./16, 3.pi./16, 5.pi./16 and 7.pi./16 etc. This method has disadvantages, however, because since the function values are not created in the desired succession a large table (of length N) must be created before the values can be output. This also involves a calculation delay before the output of the first desired values.
In a second method described in (1), the same interpolation equation is employed, but instead of creating a large table of length N holding all the desired function values, the algorithm employs a short table of length (log.sub.2 N)+2, which at any given time contains the current function value of the sequence, the two values between which the next function value of the sequence will be interpolated, and the values necessary to progressively interpolate all further function values in the sequence; whenever an entry value in that table has been used, it is replaced by a new interpolation value to be used in a later step. The nature, and the total number of calculations is the same as in the previous method, but the storage required for the table is considerably reduced and the results appear, in sequence, without initial delay.
Of course, the two methods are not necessarily mutually exclusive.
A way of implementing this latter interpolation method is as follows.
We begin by looking at the working vector or array V, of size (log.sub.2 N)+2, from which a sequence of sine values are obtained. This vector has some of the sine values which will be output as sample values. As a sample value is output, it is replaced by another sample value which will be used later on. Below is a five element vector, which shows the flow of data storage and computation in calculating up to 128 values of a sine. For simplicity, an element is shown by an integer number, which in reality represents a sine value. For example, the 16 in row 1 (at position V.sub.4) stands for sin (16*.phi..sub.o), where .phi..sub.o is the angular spacing between sine samples.
TABLE 1 ______________________________________ j = n 4 3 2 1 0 ______________________________________ 1 16 8 -4 -2 1 2 16 -8 -4 2 3 3 16 8 -4 -6 3 4 --16 -8 4 6 5 5 16 -8 12 -6 5 6 16 -8 -12 6 7 ______________________________________
Looking at the situation at the stage where the first output value is to be output (time n=1), the italicised element in that row, sin (1*.phi..sub.0), is output (V.sub.0 is output), and replaced by sin (3*.phi..sub.0). This value is computed using the underlined elements in the row, and with the interpolation formula EQU sin (K*.phi..sub.0)=[(sin ((K+L).phi..sub.0)+sin ((K-L).phi..sub.0)]/[2* cos (L.phi..sub.0)] (1a)
for K=3, L=1 in this case. In general the Ls are powers of two, and so the factors 1/[2* cos (L.phi..sub.0)] are precomputed, stored, and then looked up when necessary. The other rows, representing successive sample output times, are similarly output and computed.
A sequence of cosine values can be generated in a similar manner, since there is a cosine interpolation EQU cos (K*.phi..sub.0)=[cos ((K+L).phi..sub.o)+cos ((K-L).phi..sub.0)]/[2* cos (L.phi..sub.0)] (1b)
Another vector of the same size as V is needed for storing the cosine values, but the same j, j.sub.a, and j.sub.b apply. Thus the cosine sequence is cheap and easy to produce simultaneously with the sine sequence.
There can be an accuracy problem with both these methods, especially when the factor 1/[2* cos (L.phi..sub.0)] has a large magnitude (as will occur when L.phi..sub.0 is around .pi./2, 3.pi./2 etc). This implies that the numerator part of the interpolation, sin ([K+L].phi..sub.0)+sin ([K-L].phi..sub.0), has a small value, implying subtraction of two nearly equal quantities. The subtraction gives round-off error in finite precision hardware (more so in fixed-point than in floating point hardware). Furthermore, the inaccurate interpolation is usually itself used to produce further interpolations, so that the error propagates through further calculations.
Accordingly, in the invention, the problem can be overcome by using a different trigonometric interpolation around such argument values. According to the invention, there is provided a trigonometric function generator for generating a value of a trigonometric function in response to an input argument comprising: input means for receiving an input argument; first store means for storing values of the function corresponding to spaced argument values; and arithmetic means for interpolatively calculating, using a first algorithm, values of the function corresponding to arguments lying in the spaces between those spaced argument values; characterised in that it further comprises:
second store means for storing values of a second, complementary function corresponding to spaced argument values; PA1 arithmetic means for interpolatively calculating, using a second algorithm, values of the first function corresponding to arguments lying in the spaces between those spaced argument values from values of the second function corresponding to the neighbouring predetermined argument values; and PA1 control means for selectively employing the second store and algorithm for an input argument which would, using the first store and algorithm, give rise to a substantial round-off error in said function value.
For example, in the case of the above method, instead of computing . . . EQU sin (K*.phi..sub.o)=[(sin ([K-L].phi..sub.o)+sin ([K-L].phi..sub.o)]/[2* cos (L.phi..sub.o)] (1a) EQU (or=[(sin ([K-L].phi..sub.o)+sin ([K-L].phi..sub.o)]* sec (L.phi..sub.o)/2)
where cos (L.phi..sub.o) is very small, we use EQU sin (K*.phi..sub.o)=[(cos ([K+L].phi..sub.o)-cos ([K-L].phi..sub.o)]/[2* cos (.pi./2+L.phi..sub.o)] (2a) EQU cos (K*.phi..sub.o)=[(sin ([K-L].phi..sub.o)-sin ([K+L].phi..sub.o)]/[2* cos (.pi./2+L.phi..sub.o)] (2b) EQU and cos (.pi./2+L.phi..sub.o)=sin (L.phi..sub.o).
Hence the sine value is found by using entries from the cosine vector (and vice versa). This cooperation between the two vectors ensures that accuracy is preserved at all times. The decision between using one or other of the above formulas depends on the L.phi..sub.o value. If L.phi..sub.o is closer to n.pi. than m.pi./2, where M is odd, use equation (1a). Otherwise, use equation 2a. A numeric example will suffice to illustrate the benefit of the invention. Let K.phi..sub.o =20.degree. L.phi..sub.o =89.9.degree. and 4 place accuracy used throughout. Equation (1a) gives EQU sin (20.degree.).perspectiveto.286.5(0.9403-0.9391)=0.3438
whereas equation (2a) yields the correct result EQU sin (20.degree.)=-0.5000(-0.3404-0.3437)=0.3420
At least one decimal place of accuracy has been lost with equation (1a).
The price paid for this is the extra memory space required for storing the complementary function interpolation values and the (small number of) extra calculations per cycle needed to update these stored values. Of course, if it is in any case desired to generate the complementary function (as, for example, in Fourier analysis) then the benefit of the invention is obtained without any extra memory being required.